Laplacian-level density functionals for the 
kinetic energy density and exchange-correlation energy 
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We construct a Laplacian-level meta-generalized gradient approximation (meta-GGA) for the 
non-interacting (Kohn-Sham orbital) positive kinetic energy density r of an electronic ground state 
of density n. This meta-GGA is designed to recover the fourth-order gradient expansion r GE4 in the 
appropiate slowly- varying limit and the von Weizsacker expression t w = |Vn| 2 /(8n) in the rapidly- 
varying limit. It is constrained to satisfy the rigorous lower bound r (r) < r(r). Our meta-GGA is 
typically a strong improvement over the gradient expansion of r for atoms, spherical jellium clusters, 
jellium surfaces, the Airy gas, Hooke's atom, one-electron Gaussian density, quasi-two dimensional 
electron gas, and nonuniformly-scaled hydrogen atom. We also construct a Laplacian-level meta- 
GGA for exchange and correlation by employing our approximate r in the Tao, Perdew, Staroverov 
and Scuseria (TPSS) meta-GGA density functional. The Laplacian-level TPSS gives almost the same 
exchange-correlation enhancement factors and energies as the full TPSS, suggesting that r and V 2 n 
carry about the same information beyond that carried by n and Vn. Our kinetic energy density 
integrates to an orbital-free kinetic energy functional that is about as accurate as the fourth-order 
gradient expansion for many real densities (with noticeable improvement in molecular atomization 
energies) , but considerably more accurate for rapidly-varying ones. 
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I. INTRODUCTION 

In ground-state density functional theory, the nonin- 
teracting kinetic energy (KE) T s of a system of N elec- 
trons may be treated as an exact functional of the oc- 
cupied Kohn-Sham (KS)i orbitals {&}, and only the 
exchange-correlation energy has to be approximated. 
However finding an accurate orbital-free KE functional 2 
will simplify and speed up by orders of magnitude any 
KS self-consistent calculation. In general, a kinetic en- 
ergy density (KED) is any function which integrates to 
the non- interacting kinetic energy T s : 



T s[n^n L ] = I dr r, 



(1) 



where nj (r) and 7ij.(r) are the electron spin densities and 
r(r) = Tj([nj],r) + Tj ( [nj , r) . Because of the local ver- 
sion of the spin-scaling relation 3 - 



r CT ([n CT ],r) = (l/2)r([n = 2n CT ],r) ! 



(2) 



we will only need to show our expressions for spin- 
unpolarized systems with = rij = n/2. There are two 
important forms of the KED: one which depends on the 
Laplacian of the Kohn-Sham orbitals and follows directly 
from the Kohn-Sham self-consistent equations: 



N 

r L (r)=-(l/2)^*(r)V 2 <Mr) 

i=l 



(3) 



and another which is positive definite: 



r(r) = (1/2) 



1. 



N 

£|V<Mr)| 2 =r L (r) + ^VMr), 

i=l 



(4) 



where n(r) = Yli=i l<M r )| 2 1S the electronic density. (We 
use atomic units, with fi = m = e 2 = 1). (r L + r)/2 has 
been proposed^ as the closest analog of a classical KED. 

While a generalized gradient approximation (GGA) 
uses only n and Vn, a meta-GGA (MGGA), such as the 
accurate nonempirical TPSS 5 density functional for the 
exchange-correlation energy, is constructed from local in- 
gredients n(r), Vn(r), and r(r). In this work, we will 
present evidence that r(r) and V 2 n(r) can carry essen- 
tially the same information beyond that carried by n(r) 
and Vn(r) . To do so, and to make an improved semi-local 
density functional for T s , we shall construct a Laplacian- 
level meta-GGA for r(r). Our semi-local functional re- 
covers the fourth-order gradient expansion (GE4) KED 
for a slowly- varying density and the von Weizsacker KED 
for a rapidly- varying density. 

The gradient expansion (GE), which becomes exact for 
densities that vary slowly over spaced, is^ 



r = T^F s (p,q,...), 



with 



E^ (2n) 



(5) 



(6) 



where p = |Vn| 2 /{4(37r 2 ) 2 / 3 n 8/3 } and q = 
V 2 n/{4(37r 2 ) 2 / 3 n 5//3 } are dimensionless derivatives 
of the density, = ^(37r 2 ) 2 / 3 n 5 / 3 is the Thomas- 
Fermi KED&, and Fs 2n \p, q, ...) is the enhancement 
factor of the 2n-th term of the gradient expansion. The 
zero-th order enhancement factor is® R (0) = 1, and 



the second-order one is&i£ F, 



(2) 



(5/27)p + (20/9)?. 



The term linear in q is a key ingredient in our MGGA. 
Although this term integrates to zero in Eq. fTJ), it 



2 



iaii important for the KED. This term also indicates 
rapidly- varying density regimes (e.g., near a nucleus 
where q — > — oo and in the tail of the density where 
q — > oo). The fourth-order enhancement factor is 12 

F (i) = A = (8/81)g 2 - (l/9)pg + (8/243)p 2 > 0. (7) 

This enhancement factor is a simplified expression ob- 
tained with Green's theorem (integration by parts in Eq. 
([1])) under the assumption that n(r) and its gradients 
vanish as r — > oo. A full expression for involving 
derivatives of the density of higher than second order, is 
given in Refj 1 - . Although the integration by parts that 
leads to Eq. is inappropriate for some non-analytic 
densities 1 ^, we shall incorporate Eq. ((7} into our MGGA. 
The sixth-order terrain, even if it provides a useful correc- 
tion to the fourth-order gradient expansion for the for- 
mation energy of a monovacancy in jellium^, diverges 
for atoms after the integration of Eq. fl}, and requires 
higher derivatives of the density than we would like to 
use. For later use, we define 

F GE4 (p, q) = l+ (5/27)p + (20/% + A. (8) 

The von Weizsacker expression-- (t w = |Vn| 2 /(8n)) 
is exact for any one- or two-electron ground state, is accu- 
rate in nearly iso-orbital regions, satisfies the exact nu- 
clear cusp condition (r(0) = Z 2 n(0)/2 where Z is the 
nuclear charge}^, is exact in the r — > oo asymptotic re- 
gion (where the density matrix behaves like Eq. (8) of 
Ref.—), and gives a rigorous lower boun d 18 ' 19 

r W (r) < r(r), (9) 

or < F s . The semi- local bound of Eq. © is one of 
the most important constraints in the construction of our 
functional. The von Weizsacker enhancement factor— 

= (5/3)p > (10) 

is simple. But the von Weizsacker KE functional gives, 
in general, very poor approximate atomization kinetic 
energies (see Table III of Ref. 2 S), and this has been 
attributed 21 to its strong violation of Eq. (8) of Refill. 

Recently Tao, Perdew, Staroverov and Scuse- 
ria (TPSS)^ have constructed a nonempirical meta- 
generalized gradient approximation for the exchange- 
correlation energy. This functional, which satisfies as 
many exact constraints as a meta-GGA can (see Table 1 
of RefJ^, or Refi^ for a detailed explanation), provides 
a good overall description of atoms, molecules, solids and 
surface a 22 31 . We construct a Laplacian-level TPSS (LL- 
TPSS) by replacing r by T MGGA in the TPSS exchange- 
correlation energy per particle. 

Laplacian-level functionals for exchange and correla- 
tion have been advocated in Refi^ 2 -, and proposed in 
Ref.—. Recent interest in them is driven in part by the 
observation of Ref. 34 that V 2 n can be used along with n 
to mimic some exact exchange-correlation energy densi- 
ties. 



The paper is organized as follows. In section [TIJ we 
present the KED meta-GGA functional. In section IIII1 
we test our functional for several physical systems and 
models, and further explain its behavior. In section IW) 
we construct the LL-TPSS exchange-correlation func- 
tional, and present numerical and analytic evidence that 
V 2 n and r carry essentially the same information beyond 
that carried by n and Vn. In section |Vl we summarize 
our conclusions. 



II. CONSTRUCTION OF A META-GGA 
KINETIC ENERGY DENSITY 

Our meta-GGA for the KED is an interpolation 
between a modified gradient expansion and the von 
Weizsacker expression. There are many ways to satisfy 
the inequality of Eq. ([5]), so we will have to rely on em- 
piricism to select one of them. 

First, we construct a modified fourth-order gradient 
expansion enhancement factor 

which has the following properties: 

(1) For small p and \q\, 

F GE4-M = pGEA + Q ^ {u) 

so that it recovers the fourth-order gradient expansion 
for a slowly- varying density. 

(2) When \q\ -► oo, 

Ff Ei ~ M — » 1 + Ff , (13) 

which is the correct limit for a uniform density perturbed 
by a small-amplitude, short-wavelength density wave 3 ^. 
When p — > oo, j?ge4-m pW _|_ (jP^ which is reason- 
able for other rapid density variations. 

(3) The modified enhancement factor of Eq. (jlip sat- 
isfies a uniform damping property: 

| F GB4-M| < | F GB4|_ (W) 

This is desirable because for large p and \q\ spuriously 
large values of |f s G£;4 | can arise from truncation of the 
gradient expansion. 

We shall assume that the condition F GE4 ~ M <~ 
indicates the need for F^ IGGA = Ff '. Our meta-GGA 
interpolates between pGEi-M anc j pw The smooth in- 
terpolating function is 

( 0, z < 

/.*(*) = I ( Jj^^j f, 0<z<a (15) 
L 1) z> a : 

where < a < 1 and b > are parameters. This function 
is plotted in Fig. CD The meta-GGA KED is defined as 
follows: 

T MGGA =T^F s MGGA (p,q), (16) 



3 



where 



pMGGA _ pW _|_ IpGEi—M _ 

x/ afc (Ff^-Ff). 



F. 



w \ 



(17) 



When f ab = 1, i.e., when pGEi—M > pW + ^ Eq jjjj, 
makes F* IGGA = F GEA ~ M . When / ab = 0, i.e., when 

pGEA—M < pW^ Eq ^ makes pMGGA = pW In 

between, F^GGA j g an interpolation between F GEi ~ M 
and i 7 ^. Our positive meta-GGA KED keeps all the 
correct features of T GEi ~ M ^ tends to be exact in iso- 
orbital regions, and satisfies the important constraint of 
Eq. ©: 



pMGGA > pW 



(18) 

The meta-GGA depends on the empirical parameters 
a and b. These parameters were found numerically, by 
minimizing the following expression 

Error = i. "m.a.r.e. atoms" + 4 "m. a. r.e. clusters" 

+i "m.a.r.e.LDM(N = 8)" , (19) 



where "m.a.r.e. atoms" is the mean absolute relative er- 
ror (m.a.r.e.) of the integrated kinetic energy of 50 atoms 
and ions, "m.a.r.e. clusters" is the m.a.r.e. of the KE 
of 2e~, 8e~, 18e~, 20e~, 34e~, 40e", 58e", 92e~, and 
106e~ neutral spherical jellium clusters (with bulk pa- 
rameter r s = 3.93 which corresponds to Na) and "m.a.r.e. 
LDM(N=8)" is the m.a.r.e. of the KE of N=8 jellium 
spheres for r s — 2, 4, and 6, calculated in the liquid drop 
model (LDM) 



T. 



LDM 



(3/10)fc|iV 
a/3 



s iV 2/3 47rr? 



(20) 



where kp = (37r n) L ^ :i is the Fermi wavevector, r s 
(3/47rn) 1 / 3 is the radius of a sphere which contains on 
average one electron, n is the bulk density, and a s is the 
surface KE. The "exact" LDM value is one computed 
with the exact <j s . Because the relative errors of surface 
kinetic energies are much larger than those of the atoms 
and spherical jellium clusters, we use the LDM approach 
for calculating the jellium surface KE errors; LDM gives 
m.a.r.e. comparable to that of atoms and clusters (see 
Table []}. The densities and orbitals we use are analytic 
Hartree-Fock 3 ^ for atoms and ions, and numerical Kohn- 
Sham for clusters and surfaces (with the local density 
approximation for the exchange-correlation potential) . 

In the process of optimization, we observed that the 
m.a.r.e. decreases very slowly when the parameter b > 3 
increases, but large values of b deteriorate the KED. So, 
we chose the following set of parameters: a = 0.5389 and 
6 = 3. As we can see in Table HI this set of parameters 
gives an accuracy close to (but better than) that of the 
fourth-order gradient expansion. 

In Fig. [TJ we plot the interpolating function f a b{z) us- 
ing our choice for the parameters. Near a nucleus, there 
is a large region where Ff Ei < and thus F GEi ~ M < 0, 
making f ao — and F^ IGGA — , the correct behav- 
ior. (Nevertheless, as r — > and q — > — oo, F GE4 ~ M — > 



TABLE I: Mean absolute relative error (m.a.r.e.) of inte- 
grated kinetic energies of 50 atoms and ions, of neutral spher- 
ical jellium Na clusters (2e~, 8e~, 18e~, 20e~, 34e~, 40e~, 
58e~, 92e~, and 106e~) and of jellium surfaces (with r s — 2, 
r 3 — 4, and r s = 6) incorporated into the liquid drop model 
(LDM) for a jellium sphere with N=8 electrons (see Eq. ((20j ) . 
The atoms and ions are: H, He, Be+ 2 , Be, Be +1 , Li, Li+\ Ne, 
Ne+ 8 , Ne+ 7 , Ne+ 6 , Ar, Ar+ 16 , Ar+ 15 , Ar+ 14 , Ar+ 8 , Ar+ 6 , Zn, 
Zn+ 28 , Zn+ 20 , Zn+ 18 , Zn+ 12 , Kr, Kr+ 34 , Kr+ 26 , Kr+ 24 , Kr+ 18 , 
Xe, C+ 4 , C+ 3 , C+ 2 , N+ 5 , N+ 4 , N+ 3 , B+\ B+ 3 , B+ 2 , 0+\ 
0+ 6 , 0+ 5 , 0+ 4 , Cu, Cu+ 27 , Cu+ 26 , Cu+ 25 , Cu+ 19 , Cu+ 17 , 
Cu +11 , Cu +1 and N. 
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FIG. 1: Interpolating function f a b(z) (see Eq. (|15jl ) versus 
z, for the optimized parameters a = 0.5389 and 6 = 3. 



1.24, 



DO 

w 



l + F™, making f ab -> 1 and F* IGGA — > 1 + Ff ; 
which is at least positive and finite.) In the r 
asymptotic region, where p « q — > oo, p GEi ~ — * F & 
making f ab — * and F^ IGGA — > F™, the correct limit. 
In the slowly-varying limit, where p and q are small, 
fab(z) — > 1 , so meta-GGA recovers the fourth-order 
gradient expansion. In Fig. [21 we plot the meta-GGA 
enhancement factor versus s = ■Jp, for different values of 
q. We observe an orderly behavior of F^ IGGA for p and 
q in the range appropriate to physical densities. For the 
integrated KE, our meta-GGA is size-consistent, and sat- 
isfies the uniform scaling relation^ and the spin-scaling 
relation—. The meta-GGA nonuniform scaling behavior 
is investigated in Section [TTT1 
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FIG. 2: Enhancement factor F s MGGA (p,g) versus reduced 

gradient s — ^/p, for several values of the reduced Laplacian 

q (-10, 0, 1/3, 1, 2, 3, and 4). The parabolic asymptote is 
p w 



FIG. 4: KE integrand 4,nr t versus radial distance r for the 
Ne atom. The integral under the curve is the KE for the neon 
atom: exact: 128.546 a.u., meta-GGA: 129.312 a.u., and GE4: 
129.749 a.u.. 
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FIG. 3: KE integrand A-kt 2 t versus radial distance r for the 
He atom. The integral under the curve is the KE for the 
helium atom: exact: 2.862 a.u., meta-GGA: 2.993 a.u., and 
GE4: 2.963 a.u.. 
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FIG. 5: KE integrand 4nr r versus radial distance r for 
the 10-electron (hydrogenic orbitals) ion with nuclear charge 
Z=92. The curves EXACT S and MGGA S show the contribu- 
tion of the four s-electrons (see Ref.— ). 



III. RESULTS: KINETIC ENERGY AND ITS 
DENSITY 

In Figs. [3] and 0] we plot the integrand of the KE, 
47rr 2 T, versus radial distance from the nucleus, for the 
He and Ne atoms. The KED of the fourth-order gradient 
expansion, t GE4 , is negative near the nucleus and is not a 
good approximation for r. The meta-GGA KE integrand 
is much improved near the nucleus, and everywhere fol- 
lows very nicely the exact behavior. 

Let us consider an ion model with 10 electrons which 
occupy the first hydrogenic orbitals and with a nuclear 
charge Z—92. In such a closed-shell system, the KED is 
determined by the density of the s-electrons alon o 38 ' 39 . 
In Fig. 03 we show the exact and the meta-GGA kinetic 
energy densities for this system (for the s-electrons, and 



for the whole system, separately). 

The KED of the one-electron Gaussian density is 
shown in Fig. [S] This system does not have a cusp near 
the center (r = 0, where p — and q ~ —0.5) and, in this 
sense, it is an important test for our meta-GGA, which 
we find to be as accurate here as it is for the He atom. 

The Hooke's atom is a simple system of two interacting 
electrons in a harmonic potential. For this system, the 
exact KED is the von Wcizsacker one. The correlated 
wavefunction and its density can be calculated exactly^ 
for special values of the spring constant k. The low- 
correlation case^i corresponds to k = 0.25a.u., and the 
high-correlation case^ to k = 3.6 x 10~ 6 a.u.. A modeled 
density very similar to the exact Hooke's atom density 
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r 



FIG. 6: KE integrand 47rr 2 r versus radial distance r for the 
one-electron Gaussian density. The integral under the curve 
is the kinetic energy: exact = 0.750 a.u., meta-GGA = 0.778 
a.u., and GE4 = 0.865 a.u. 
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FIG. 7: KE integrand 4-7rr 2 r versus radial distance r for the 
pseudo-Hooke's atom in the high-correlation case, using the 
scaled density given in Eq. (|22|l . The integral under the curve 
is the kinetic energy: exact = 1.115 a.u., meta-GGA = 1.264 
a.u., and GE4 = 1.185 a.u. 



(with k = 3.6 x 1CT 6 a.u.) is 

n{r) = A(l + B C r 2 )e- Cr2 , (21) 

where A = 0.67x 1(T 6 a.u., B=ll a.u., and C=0.0010212 
a.u. Applying uniform scaling n(r) — ► 7 3 n(7r) = n 7 (r), 
we define the following sealed density: 

n 7 (r) = 0.02145(1 + 10.5r 2 ) e - r2 , (22) 

with 7 = 31.753. In Fig. \7\we show the exact, meta-GGA 
and fourth-order gradient expansion KE integrands for 
the pseudo-Hooke's atom (using the scaled density given 
in Eq. ([22]) ) in the high-correlation case. The region 
near the nucleus of the pseudo-Hooke's atom is unusual 
because there strong correlation creates a deep "hole" in 
the density: q decreases smoothly with increasing r (from 
19.29 at r = to sw at r = 0.68), and p increases slowly 
with increasing r (from at r — to a peak value of 
1.829 at r = 0.233). This region can not be described 
accurately by the meta-GGA, as we can see in Fig. [7] 
Near r = 1, where p ~ and q ~ —0.495, the KED 
meta-GGA recovers the exact behavior. 

Figs. [7|and[8]show how our T MGGA can sometimes fail 
to recognize iso-orbital regions where r = r w . For the 
2e~ jellium cluster between r = 0.1 and r = 1, where p in- 
creases very slowly with increasing r (from p = 0.0016 at 
r = 0.1 to p = 0.015 at r = 1) and q ~ -0.29, the meta- 
GGA functional switches from the fourth-order gradient 
expansion to the exact behavior due to the construction 
of the interpolating function f a b. This feature is impor- 
tant for molecules, because a similar case arises at the 
center of a diatomic molecule X^. 

In Fig. [9] we plot the KED r versus the scaled distance 
£ for the Airy gas, which is the simplest model of an edge 
electron gas^. In this model, the non-interacting elec- 
trons see a linear potential. In the tail of the density, 
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FIG. 8: KE integrand 4ivr 2 T versus radial distance r for the 
2e~ Na jellium spherical cluster. The integral under the curve 
is the kinetic energy: exact = 0.114 a.u., meta-GGA = 0.120 
a.u., and GE4 = 0.124 a.u. 



the KED meta-GGA becomes exact (see Fig. [9]) as we 
mentioned in Sec. [TTJ At least two orbital-free kinetic 
energy functionals have been based upon the Airy gas 
model. In Ref. 44 , the kinetic energy density of the Airy 
gas is transferred to other systems in a Local Airy Gas 
(LAG) approximation, which seems accurate for jellium 
surfaces but makes r diverge at nuclei. In Ref4^, a den- 
sity functional is constructed for the linear potential; it 
has an unphysically rapid oscillation 46 in its correction 
to t(°) for a slowly-varying density. 

For a quasi-2D electron gas (quantum well) whose or- 
bitals are those of noninteracting electrons confined by 
infinite barriers 4 ^, the exact KE per electron can be cal- 
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FIG. 9: Kinetic energy density versus Airy scaled distance 44 
for the Airy gas. The exact KED is given by Eq. (11) of 
Ref.-**. 
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FIG. 11: Integrated KE versus scaling parameter A for the 
nonuniformly-scaled hydrogen atom, in the prolate case. 
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FIG. 10: Integrated kinetic energy per electron versus L/rl D 
for a quasi-2D electron gas in the infinite barrier model. The 
exact curve is given by Eq. (|23[) . rl° is the areal density 
parameter of the quasi-2D electron gas and L is the width of 
the quantum well. (L < 3.85r? D , see Ref4*) Note that rj 4) 
diverges here. 



culated analytically, and has the simple expression: 



N 



IT 

2L 2 



„2D\2 
"F J 



(23) 



2D = V2/r 2 s D 

„1D 



where L is the width of the quantum well, k F 
is the two-dimensional Fermi wavevector and r 2U is the 
radius of the circle that contains on average one electron 
of the quasi two-dimensional gas. Fig. [TO] is for the elec- 
tron gas at r 2D = 4. diverges due to the infinite 
barriers, so it is not shown in the figure. As we can see, 
the KED meta-GGA functional solves this nonuniform 
scaling problem almost exactly. 



The nonuniform scaling inequality** is 
T s [nl}<\ 2 T*[n}+Ty{n}+T* 
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FIG. 12: Integrated KE versus scaling parameter A for the 
nonuniformly-scaled hydrogen atom, in the oblate case. 



where A is a positive scale factor, the nonuniformly- 
scaled density is rv^{x,y,z) = \n{Xx,y,z), and T^[n] = 
(l/2)(dT s [n q x ]/dX)\ = i, where q is x, y or z. For the 
von Weizsacker functional, and thus for an exact treat- 
ment of the nonuniformly-scaled hydrogen atom of den- 
sity n x (r) = (A/ 7 r) exp(-2 v /(Ax) 2 + y 2 + z 2 ) , Eq. (JH 
becomes an equality**. In Figs. [TT] and [TSJ, we show 
the KE of the nonuniformly-scaled hydrogen atom as a 
function of A, for the prolate case (A < 1) and oblate 
case (A > 1) respectively. As we can see, the meta-GGA 
KE functional does not satisfy the nonuniform-scaling 
inequality, but is still very close to the exact behavior. 
The von Weizsacker KE is always less than or equal to 
the meta-GGA functional (see Eq. (fT8| ). So, the meta- 
GGA seems to describe the nonuniform scaling relation 
with considerable fidelity, a potentially important fea- 
ture for molecules where the bonding causes nonuniform 
distortions in the density. 

When a molecule at its equilibrium geometry is bro- 
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ken up into separate atoms, the total energy increases 
and (as suggested by the virial theorem) the kinetic en- 
ergy decreases. In Table |n] we present the atomization 
kinetic energies for a set of molecules used in Ref^S. The 
meta-GGA kinetic energy functional gives the best over- 
all results, but is still not accurate enough for chemical 
applications. We observe that for NO and O2 molecules, 
where the other listed functionals fail badly, the meta- 
GGA keeps the right sign. (In fact, with the exception of 
N 2 and F 2 , the atomization kinetic energies of the meta- 
GGA all have the right sign.) 



TABLE II: Integrated atomization kinetic energy ( KE atoms 
- KE molecule, in a. u.) for a few small molecules. The 
kinetic energies were calculated using the PROAIMV code 
with Kohn-Sham orbitals given by the Gaussian 2000 code 
(with the uncontracted 6 — 311 + G(3df, 2p) basis set, Becke 
1988 exchange functional, and Perdew-Wang correlation 
functional 50 ). The last line shows the mean absolute errors 
(m.a.e.). 





T->cxact 
L s 


T (0) 
- 1 S 






rpMGGA 


H 2 


-0.150 


-0.097 


-0.114 


-0.119 


-0.216 


HF 


-0.185 


-0.305 


-0.186 


-0.133 


-0.352 


H 2 


-0.304 


-0.308 


-0.136 


-0.057 


-0.634 


CH 4 


-0.601 


-0.737 


-0.571 


-0.498 


-1.036 


NH 3 


-0.397 


-0.231 


-0.060 


0.014 


-0.477 


CO 


-0.298 


-0.323 


-0.085 


0.015 


-0.458 


F 2 


-0.053 


0.128 


0.282 


0.338 


0.154 


HCN 


-0.340 


-0.1835 


0.079 


0.186 


-0.328 


N 2 


-0.158 


0.344 


0.565 


0.650 


0.319 


CN 


-0.431 


-0.215 


0.005 


0.094 


-0.231 


NO 


-0.268 


0.092 


0.330 


0.422 


-0.084 


O2 


-0.100 


0.106 


0.335 


0.431 


-0.194 


m.a.e. 




0.177 


0.311 


0.384 


0.201 



IV. LAPLACIAN-LEVEL META-GGA FOR 
EXCHANGE-CORRELATION ENERGY 

The exchange-correlation meta-GGA uses as ingredi- 
ents the spin densities and their gradients, and the posi- 
tive Kohn-Sham kinetic energy densities: 



xF, 



MGGA 



E mgga = Jdrn(r)er{n(r)) 
(nr (r), ni (r), Vn T (r), Vn T (r), r T (r), n (r$5) 



where e£\n) = (— 3/ An) (Sir 2 n) 1 ^ 3 is the exchange energy 
per electron of a uniform, spin-unpolarized density n. 
The local density approximation (LDA) is recovered for 



Vri| = Vri| = and t\ = T \°\ T i = T \ V '' ■ Our Laplacian- 
level non- interacting KED functionals r^ IGGA (r) and 
r^ IGGA {r) can replace the exact Kohn-Sham kinetic en- 
ergy densities r-f (r) and r^(r)) in Eq. (|2"5")) . The resulting 
LL-MGGA is an explicit density functional, while Eq. 
(I25P is an implicit one. 



.(0) 




FIG. 13: Exchange enhancement factor F x versus 
reduced gradient s = ^/p, for several values of the reduced 
Laplacian: q = —10, —1/3, 0, 1/3, 1, 2, 3. 



We have constructed in this way the Laplacian-level 
(LL-TPSS) of the Tao, Perdew, Staroverov and Scuseria 
exchange-correlation meta-GGA 5 .. Like TPSS, LL-TPSS 
recovers the fourth-order gradient expansion 5 ^ 5 ^ for the 
exchange energy in the slowly-varying limit. Moreover, 
like TPSS, LL-TPSS has a finite exchange potential at 
a nucleus. The need for an ingredient beyond n and V?i 
(e.g., V 2 n) to satisfy this exact constraint was empha- 
sized in Refsi22 and^. The nuclear cusp of an atom can 
be defined by q — > —00 and s ~ 0.376, so the constraint 
used in the construction of the TPSS exchange enhance- 
ment factor 5 



dF^ PSS (s,Z^ l)/ds| s=0 .376 = 0, 

where z = t(y)/t w (y), becomes 

dF^ L - TPSS (s,q^ -oo)/d S | s=0 . 376 = 0. 



(26) 



(27) 



Such constraints can be satisfied by a Laplacian-level 
meta-GGA, but not by a GGA (using only n and Vn). 

In Fig. we show the exchange enhancement factor 
F^ L ^ TPSS as a function of the inhomogeneity parame- 
ter s = y/p for several values of the reduced Laplacian 
q. The enhancement factor interpolates in an orderly 
way between the exact slowly- varying limit (for small s 
and |<7|) and a rapidly-varying limit (for large s), while 
satisfying Eq. (|27|) . 

In Figs. Hi] and US] we compare the TPSS and the LL- 
TPSS enhancement factors (for exchange and exchange- 
correlation) for the Zn atom, which was also studied in 
Ref.— . The LL-TPSS exchange-correlation energies are 
close to the TPSS values for this and other atoms. (For 
the H atom, E^ PSS = and E LL - TPSS = -1.537 -lO" 6 ; 
for the other 49 atoms and ions, the m.a.r.e. of LL-TPSS 
with respect to TPSS is 0.00132 for exchange, 0.0098 
for correlation, and 0.00133 for the combined exchange- 
correlation energy.) 
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FIG. 14: Exchange enhancement factors F^ L TPSS and 
pTPSS versus rac ii a i distance r for the Zn atom. E^ PSS — 
-69.798 a. u., E^ L ~ TPSS = -69.528 a.u.. Note that 
< r _1 > _1 is 0.65 for the 3d and 2.26 for the 4s electrons. 



FIG. 16: Deviation (from local density approximation) of the 
exchange energy integrand 4irr 2 e x versus radial distance r for 
the Xe atom. GEA is the second-order gradient expansion 
for exchange. E^ PSS = -178.449 a.u. and E LL - TFSS = 
-178.240 a.u.. 
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FIG. 15: Exchange-correlation enhancement factors 
Fxc~ TPSS , F^ C PSS , F P C BE and F E C DA versus radial distance 
r for the Zn atom. E^ PSS = -71.208 a.u., E LB ~ TPSS = 
-71.073 a,.u.,E P c BE = -70.934 a.u.. Note that < > _1 
is 0.65 for the 3d and 2.26 for the 4s electrons. PBE is the 
non-empirical GGA of Ref.— 

In Fig. [THwe show that the LL-TPSS exchange func- 
tional, like the TPSS, shows a strong enhancement in 
the Is region of an atom but is elsewhere not so different 
from the second-order gradient expansion for exchange, 
as discussed in Ref£. 



V. CONCLUSIONS 

We have constructed a Laplacian-level kinetic energy 
density (KED) meta-GGA functional which depends on 
two empirical parameters. These parameters control an 
interpolation between a modified fourth-order gradient 
expansion and the von Weizsacker expression. This in- 



terpolation is designed to respect the exact constraint of 
Eq. ©. From our tests and results of section IfTTl we con- 
clude that the KED meta-GGA functional is a successful 
approximation to the positive Kohn-Sham KED. 

Our functional uses a simplified expression for the 
fourth-order gradient expansion^ of r. We also tested 
the meta-GGA with the full fourth-order terms^, but 
did not find a considerable improvement in KED or in- 
tegrated KE. 

We have built a Laplacian-level meta-GGA for ex- 
change and correlation (LL-TPSS) which seems to im- 
itate faithfully the TPSS^ exchange-correlation meta- 
GGA. Several constraints of the TPSS meta-GGA are 
exactly satisfied by the LL-TPSS, but others, like E c = 
for any one-electron system, are approximately satisfied. 
It appears to us that V 2 n and r carry essentially the 
same information beyond that carried by n and Vn. Ei- 
ther V 2 n or t can be and are used to recover the fourth- 
order gradient expansion of the exchange energy in the 
slowly-varying limit and to make the exchange potential 
finite at a nucleus. 

Integration of our Laplacian-level meta-GGA for the 
kinetic energy density yields an orbital-free density func- 
tional for the kinetic energy that seems to improve 
upon the fourth-order gradient expansion, especially for 
rapidly- varying densities. We have made all our tests for 
electron densities constructed from orbitals, and do not 
know what might be found from a selfconsistent solution 
of the orbital- free Euler equation for the electron density. 
It was argued in Ref^ that KE functionals employing 
only n and Vn cannot yield both accurate integrated en- 
ergies and accurate functional derivatives. We believe 
that accurate results could be expected in some rapidly- 
varying regions (near nuclei and in density tails) from KE 
functionals that employ n, Vn, and V 2 n; our specific ex- 
pressions however may encounter a problem due to the 
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sharper features in Fig. [2] (also visible in Fig. [T3|) . It 
is generally believed that correct quantum density oscil- 
lations and shell-structure oscillations can only be found 
from a fully-nonlocal orbital-free density functional for 
the kinetic energ y 54 i 55 . 
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